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slowly  when  confronted  by  the  large  representations  necessary  for  low-level 
vision.  Application  of  muligrid  methods  can  overcome  this  drawback.,  as  we  establish 
in  previous  work  on  3-D  surface  reconstruction.  In  this  paper,  we  develope 
efficient  multiresolution  iterative  algoithms  for  computing  lightness,  shape- 
f rom-shading,  and  optical  flow,  and  we  evaluate  the  performance  of  these  algorithms 
on  synthetic  images.  The  multigrid  methodology  that  we  describe  is  broadly 
applicable  in  low-level  vision.  Notably,  it  is  an  appealing  stategy  to  use  in 
conjuntion  with  regularization  analysis  for  the  efficient  solution  of  a 
wide  range  of  ill-posed  visual  reconstruction  problems. 


,  D  ■  sir i,  utli  i\,‘ 

I  Av*;llftbl  1  :  '  V 
\-  -  - 

i  (Avail  t* /  •: v 

Spoc  :  ij. 


|Dl3t 


MASSACHUSETTS  INSTITUTE  OE  TECHNOLOGY 
ARIII ICiAl  INI  El  LICENCE  LABORATORY 


A.].  Memo  No.  803 


October,  1584 


THE 


MULTIGRID  RELAXATION  METHODS  AND 
ANALYSIS  OF  LIGHTNESS,  SHADING,  AND  FLOW 


Demetri  Terzopoulos 


Abstract 

Image  analysis  problems,  posed  mathematically  as  variatir  ml  principles  or  ns  partial  differential 
equations,  are  amenable  to  muneiieal  solution  by  relaxation  algorithms  that  are  local,  iterative, 
ar.d  often  parallel.  Although  they  arc  well  suited  structurally  for  implementation  on  massively 
parallel,  !i  eally  inteiumiiecud  computational  architectures,  such  distributed  algorithms  are  seriously 
handicapped  by  an  inherent  inefficiency  it  propagating  constraints  between  ^idely  separated 
processing  elements,  lienee,  tbey  converge  extremely  slowly  when  confionlcd  by  the  large 
representations  necessary  for  low-level  vision.  Application  of  multigrid  methods  can  overcome 
fa:s  drawback,  as  we  established  in  previous  work  on  .3-1)  surface  reconstruction.  In  this  paper,  we 
develop  eflicicnt  ir.ultoesolution  iterative  algorithms  for  computing  lightness,  shape-from-shading, 
and  optical  lltjv. .  and  we  evaluate  die  performance  of  these  algorithms  on  synthetic  images.  The 
u.uitig’.i iJ  mctln,  .ology  Out  we  describe  is  broadly  applicable  in  low-level  vision.  Notably,  it  is  an 
appealing  strategy  to  use  in  conjunction  with  regularization  analysis  for  die  eflicicnt  solution  of  a 
wide  range  of  ill-posed  ’.nun!  reconstruction  problems.  ,  ' 


;'c;  Massachusetts  Institute  of  Technology  1984 


l  iu-  i <  .tevribp<  i  •'e-cch  done  at  the  Aitificial  Intelligence  Laboratory  of  the  Massachusetts 
Institute  i  i  i  cchnul  >gy.  Support  for  the  laboratory’s  Artificial  intelligence  research  is  provided 
in  pint  In  the  Advanced  Kese.ucli  Projects  Agency  of  lire  Department  of  Defense  under  Office 
O'  N.o.'l  Tiscaich  comi.m  NtJi’di  t-yiJ  C-0505  and  the  System  Development  Foundation.  The 
nuihoi  w.a:  supported  by  -he  Natui.d  S<  wnr  and  Engineering  Research  Council  of  Canada  and 
the  l  onds  I  d’.A.C  .  OuibN,  Canada 


1.  Introduction 


Variational  principles  and  partial  differential  equations  have  played  a  significant  role  in 
the  mniheir, alien!  formulation  of  low-level  usual  information  puKessing  problems  (reprcsctuaiivc 
examples  include  [Horn,  1974,  1975;  Gilman,  1979;  Horn  &  Schunck,  1981;  Ikeuchi  &  Horn; 
1981;  Narayanan  et  tii..  1982;  llajesy  &  Broil.  1982;  Hummel  &  /.ucker.  1983;  Grimson.  1983; 
Tcrzopoulos,  1982,  1983;  Nagel,  1983;  Hildreth,  1984;  Brady  &  Yuillc,  198-1JV  attractive  feat 
of  vuriatii  nal  and  differentia!  formulations  (once  discretized)  is  die  possibility  of  computing  l  c 
desir’d  solutions  by  a  popular  class  of  numerical  relaxation  algorithms.  These  iterative  algorithms 
require  oily  local  computations  which  can  usually  be  performed  in  parallel  by  many  locally 
communicating  processors  distnhuieJ  in  computational  networks  or  grids. 

I.ocil.  parallel  algorithms  arc  appealing  in  die  context  of  lovv-icvcl  \ ision  [Roscnfeld  cl 
1976;  Ldlman.  1979;  Ballard  ci  ill.,  1983J.  At  a  ceitain  level  of  abstraction  dicy  do  not 
appear  in*  ompaiiblc  with  die  apparent  structure  of  advanced  biological  vision  systems.  Moreover, 
they  are  iicaily  suited  to  implementation  on  massively  parallel  computers  with  numerous  simple, 
locally  interconnected  processing  elements.  Such  potentially  powerful  architectures  will  certainly 
proliferate,  pending  imminent  advances  in  VLSI  technology  [Batcher,  1980;  Ilillis,  1981]. 

flic  desired  solutions  to  many  visual  problems  appear  to  possess  certain  global  properties 
%  •  (consistency,  smoothness,  minimal  energy,  etc.),  winch  arc  expressed  formally  by  die  variational 

principle  or  associated  pmtial  differential  equation  formulations.1  Given  only  local  communication 
capabilities  among  processing  elements,  however,  global  properties  can  only  lie  satisfied  indirectly , 
typically  by  iteratively  propagating  visual  constraints  across  the  grid  network.  Indirect  propagation 
can  result  in  substantial  computational  inefficiency,  since  die  computational  grids  necessary  for  low- 
level  vision  applications  tend  to  be  extremely  huge.  Convergence  of  die  iterative  process  is  often 
sn  slow  as  to  nearly  ncutralt/e  die  computational  power  offered  by  massive  parallelism.  Indeed, 
foi  fine  discreti/atioiis  on  huge  grids,  excruciatingly  slow  convergence  rates  have  been  observed  in 
iterative  algmithms  for  computing  lightness  [Blake,  1984,  sec  also  Horn,  1974],  shape-from-ahading 
[Ikeuchi  cY  Horn.  1981.  Smith  I9ji?|.  optical  Haw  [ 1 1  orn  k  Scliunek.  1981;  Nagel,  198)].  3-1) 
surfaces  [Giimson,  1983;  IVtzopmilos.  1982,  1983).  and  other  visual  reconstruction  problems. 

Sme:  spatial  locality  of  computation  is  dependent  on  spatial  resolution,  local  (c.g.,  nearest 

iieieliboi)  compulations  on  a  con: sc  grid  o\ci  a  given  region  arc  analogous  to  more  global 

compuiaimis  on  a  fine  giid  over  the  same  region.  This  suggests  die  possibility  of  counteracting 

the  sluggishness  of  global  mtei actions  by  deploying  local  iterative  processes  over  a  imiltiiesolution 

hiciaichy  of  ends.  This  is  die  basis  of  million. I  rehiciitiim  methods  which  are  gaining  popularity 

in  applied  mime:  icul  analysis  |  Hack  bust  h  &  i  rottciherg,  1982).  The  computational  structure  of 

’  Variation  il  and  dill'crcminl  ibminhtUoiis  can  be  relakil  llimuuli  the  f  uler  I  agrnnee  equations  of  die  calculus 
ut  va.iatioi  s.  Eivc-n  appr.ipoaie  co:i!in:::iv  and  svnimetri  (or  self  adjoint ''ess)  conditions  [C'our.int  &  llilbeil, 
1953], 
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mulligrid  methods  bears  an  mtcicstmg  analogy  lo  the  mulliicsoltiuon  namic  of  spatial  frequency 
channels  in  tiie  human  early  visual  system  [Bruddick.  cl  a!..  1978].  The  methods  arc  also  related 
to  certain  mulliresolution  image  processing  structures  lliat  have  been  proposed,  notably  pyramids 
(Rosen  fold,  1984], 

In  earlier  work,  we  developed  an  efficient  smfacc  reconstruction  algorithm  based  on  mulligrid 
relaxation  methods  [Ter/opoulox.  1982.  l‘>33]  and  we  suggested,  as  has  Cila/cr  [1984],  that  mulligrid 
methods  are  broadly  applicable  in  low-ievc!  computer  vision.  After  a  brief  overview  of  multigrid 
methodology,  we  apply  it  to  three  other  vision  problems:  the  well-known  problems  of  computing 
lightness,  shape -from-shading.  and  optical  flow  from  images.  We  develop  novel  multircsolution 
algorithms  for  each  problem.  On.  empirical  results  indicate  that  these  algorithms  offer  ordcr-of- 
macnitudc  gains  in  efficiency  over  their  conventional  single  level  counterparts. 


2.  Multigrid  Methodology 

Pkmecimg  investigations  into  multigrid  methodology  include  the  work  of  (Fedorenko,  1961], 
(Bakhvalov,  1966],  (Brandi,  1773,  1977],  and  (Nicolaidcs.  1977],  It  has  been  applied  to  many 
boundary  vjlwe  problems  (see  (Brand,  1982]  for  an  extensive  bibliography)  and  tlicre  has  also  been 
some  development  in  the  context  of  variational  piobiem  [Nicolaidcs,  1977;  Biandt,  1980], 

2.1.  Multigrid  Relaxation  Methods 

Multigrid  relaxation  methods  take  advantage  of  multiple  discretizations  of  a  continuous 
problem  over  a  range  of  resolution  levels.  The  comscr  levels  trade  off  spatial  resolution  for 
direct  communication  paths  over  larger  distances.  Hence,  they  effectively  accelerate  the  global 
propagation  of  information  to  amplify  die  overall  efficiency  of  the  iterative  iclaxation  process. 

The  inherent  computational  sluggishness  of  local  iterative  algorithms  can  be  studied  from  a 
spatial  frequency  perspective.  A  local  homier  analysis  of  die  error  function  (or,  more  conveniently, 
tiie  dynamic  residual  function)  from  one  delation  to  die  next  shows  that  high-frequency  components 
of  the  error  —  those  components  with  wavelengths  on  the  order  of  the  grid  spacing  —  arc  short¬ 
lived.  wheieas  low- frequency  components  persist  through  many  iterations  (Brandi  1977],  Hence, 
common  or  cron  norms  dec i ease  sharply  during  die  first  few  iterations,  so  long  as  there  are 
high-frequency  components  to  be  annihilated,  but  soon  degenerate  to  a  slow,  asymptotic  diminution 
when  only  low- frequency  components  remain  (see  Fig.  1).  Ill  is  suggests  Unit  while  relaxation  is 
incllicicm  at  completely  annihilating  the  enor  function,  it  can  be  very  ellicicnt  at  smoothing  it. 

!  rom  this  point  of  view,  the  end  hierarchy  enables  die  etliuent  smooihi.v  pn  miles  of  relaxation 
to  he  exploited  over  a  wide  range  of  spatial  frequencies. 

empirical  studies  of  model  problems  (Poisson's  equation  in  a  rectangle)  indicate  dial  multigrid 
methods  can  converge  in  essentially  order  (>(  V)  number  i  t  operations,  whole  A  is  the  number  of 
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Figure  1.  Asymptotic  error  reduction  by  a  nx.uion  Hie  mean  square  (dynamic  residual)  error  is  plotted 
.is  a  tuna  on  of  the  iteration  number  tor  a  sequence  of  Kj.tuss-Seidcl)  relaxation  iterations  of  a  surface 
recoitstrueiion  algorithm.  Cue  curve  exhibits  a  typical  behavior  of  local  iuiativc  methods:  Convergence  is 
rapid  during  the  first  lev.  iterations,  but  quickly  degenerates  to  slow  asymptotic  error  reduction. 

nodes  in  the  g.'  id  [Hr, mat.  !97"|.  This  car,  be  compiled  to  typical  complexities  of  0(.VJ)  operations 
for  the  solution  of  mode!  problems  by  standard  (single  level)  relaxation.  As  a  consequence, 
•  nudtierid  methods  potentially  offe:  dramatic  increases  in  efficiency  over  standard  relaxation  methods 

in  low -lex el  vision  applications,  since  A'  tends  to  be  very  large  (order  1  o 1  to  10s,  or  more).  Tor 
comparative  complexity  analyses,  die  total  computational  expense  of  multigrid  methods  may  be 
measured  in  convenient  machine  independent  units,  lire  basic  work  unit  is  defined  ns  die  amount 
of  computation  requited  to  perform  one  iteration  on  the  finest  grid  in  the  hierarchy. 

Ot u  adaptation  of  uiultigrij  methods  to  visual  processing  has  a  number  of  features:  (i) 
multiple  visual  represent. -.dons  covering  a  range  of  spatial  resolutions,  (ii)  local,  iterative  relaxation 
processes  dr.'ll  propagate  constraints  within  each  representational  level,  (lit)  local  coarsc-to-finc 
I'raiony.iition  processes  that  allow  coarser  representations  to  constrain  finer  ones,  (is)  frnc-to-coarsc 
irs!r;ii:ni:  processes  that  allow  finer  representations  to  constrain  and  improve  the  accuracy  of  coarser 
('ties,  uul  (tv)  (iccuisivc)  eooidin.iiian  schemes  that  enable  the  liioi. uchy  of  representations  and 
component  pro  cesses  u,-  cooperate  h  wards  increasing  cllieieney. 

In  r  itiltigrid  methods,  the  mti.ikvel  processes  usually  arc  basic  relaxation  mediods  such  as 
Ciauss-Sc-.dcl  or  Jacobi  relaxation,  die  prolongation  prtKesscs  arc  local  Lagrange  (polynomial) 
iiiteipol.U'ons,  and  tire  rest r icimr,  processes  are  local  averaging  operations.  1  lie  exact  fomi  of  these 
c'peratioii ,  is  problem  depeiulcm. 


2.2.  Discretization 

Appropriate  relaxation  processes  can  be  derived  by  local  discretization  of  die  continuous 
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problems.  I ~hc  finite  dement  method  |Str.mg  &  fix,  1973],  a  general  and  powerful  local 
discretization  technique,  can  be  applied  directly  to  variational  principle  foi initiations  of  visual 
problems  [1  cr/opoulos.  1982],  When  the  visual  problem  is  posed  as  a  partial  differential  equation, 
lucul  discretization  may  be  carried  out  using  the  finite  difference  method  [l  orsy the  &  Wasow,  I960], 

The  basic  idea  behind  the  finite  element  method  is  that  a  global  approximation  can  result  from 
interactions  among  many  very  simple  local  approximations.  This  is  accomplished  by  tcssellating 
the  continuous  domain  into  a  large  number  of  small  subdomains  or  elements  If  whose  dimensions 
depend  on  a  fundamental  size  h.  I  he  appioxnn.iiion  within  elements  depends  on  a  small  number  of 
parameters  —  the  values  of  the  solution,  and/or  some  of  us  derivatives,  at  a  set  of  nodes  associated 
with  each  clement.  The  power  of  the  method  stems  from  the  fact  die  local  approximations  can  be 
based  on  low-ordci  polynomials.  Ibis  makes  it  relatively  easy  to  express  the  continuous  functional 
as  a  discrete  summation  over  all  the  element  contributions.  If  the  variational  principle  is  quadratic, 
the  resulting  discrete  problem  takes  the  form  of  a  large  system  of  linear  equations  Ahu*  =  f\ 
whcic  u*  is  the  vector  of  nodal  variables.  The  finite  element  method  can  also  he  characterized  as  a 
systematic  procedure  for  generating  finite  element  approximating  spaces  whose  iocal-support  basis 
functions  make  V*  sparse  (i.c.,  most  of  its  elements  arc  zero). 

The  finite  difference  method  is  applied  differently.  Typically  a  giid  of  nodes  with  spacings 
proportional  to  a  parameter  h  is  set  up  over  the  domain.  1  he  diltcreiui.il  operator  is  then  replaced 
by  finite  dilfercnco  equations  involving  nodal  variables  at  neighboring  nudes.  The  collection 
of  finite  difference  equations  defines  a  disci  etc  system  which  approximates  the  given  differential 
equation.  If  die  differential  operator  is  linear  (as  aie  the  Hulcr-I  agrange  equations  of  quadratic 
variational  principles)  and  a  linear  finite  difference  approximation  is  employed,  the  discrete  system 
is  again  a  linear  system  Ahuh  —  f\  Although  the  total  number  of  nodes  ,V  is  generally  large,  each 
finite  difference  equation  involves  only  a  few  nodal  variables.  Theieforc.  the  linear  system  is  again 
sparse. 

While  the  finite  difference  method  is  generally  easier  to  apply,  the  finite  clenwni  method  offers 
a  much  sounder  convergence  theory,  as  well  as  a  flexibility  tliat  allows  the  spatially  nonui.iform 
discretization  of  domains  having  complicated  shapes.  Nonetheless,  both,  discieti/ation  techniques 
yield  large,  sparse  systems  of  linear  equations  in  a  wide  range  of  visual  applications.  A  great  deal 
of  ell’ort  in  nume:  wai  .nialvsix  has  been  directed  to  the  solution  of  such  systems,  which  turn  out  to 
be  especially  well  suited  for  solution  by  local,  p.u.dlcl.  iterative  methods,  particularly  die  relaxation 
methods  that  wc  have  been  discussing. 

2.3.  Multigrid  Struclure  and  Coordination 

Our  spatially  uniform  discretizations  of  the  continuous  vi-ual  problems  in  this  paper  will 
yield  uniform  grids  at  e.,Ji  level  of  die  mulligiid  hie:  ueby.  Application  of  mulligiid  methods 
can  be  simplified  substantial!,  giver,  a  2:1  dec  lease  in  god  lesuh.-tion  fio:n  any  level  to  the  next 
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figure  2.  PomHe  gini  oie.iiii/unnn  of  ,i  nv.. a. resolution  algorithm.  A  small  p' >r lion  of  three  levels  of  the  2:1 
inultiei kI  hiemrciiv  t>  Or.h  ncdiovnciglihnr  in’eiprocessor  connections  .ire  included. 

coarser  Ic'.cl.  l  ot t'lo.iu !y.  tins  icvifition  ratio  appears  to  be  near  optimal  with  regard  to  multigrid 
coincidence  rates  (lii.oult,  197~|.  lag.  2  illustraics  a  portion  of  three  grids  of  a  2:1  multigrid 
’Hierarchy.  In  a  serial  implement. iimn  the  eenual  pmcessor  operates  at  each  grid  node  sequentially, 
whereas  in  a  fulls  pan  di.-i  m.nlemmiucon.  each  node  represents  a  separate  processing  element 
vilhin  a  distributed  '.oml-mteicuiUKc!  e.Khiisctiuc  (see  l  ig.  2). 

I  he  lmiltiicsoluaon  vim;.?!  ile.uithms  to  no  desenhed  utili/c  sim pic  injection  for 

the  line- tt -coarse  resuimions.  bilinear  liuerpolaiion  for  the  coarsc-to-finc  prolongation, 

and  .in  tii.-.ij.iivc  multigrid  uordmation  scheme  which  was  employed  successfully  in  our  surface 
reconstitution  afmntlui  (see  J I  er/opoulos.  WS2.  1 983]  for  details).  flic  general  coordination 
1  me  lira  performs  a  v.illi-.icm  number  of  relax. nion  iterations  to  solve  the  coarsest  level  discrete 
Sjj.cm  A h  ii'*1  P-  in  desired  .u curacy  (procedure  SOLVE),  and  then  priced*  to  the  finest  level 
l  I ,  atet  i ding  tt> 
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procedure  FMG 

uKl  -  SOLVE  (1,  uh‘,  fh| ) ; 
for  1^2  to  A  do 
beg  i  n 

>h‘  ♦-  Ii-i ; 

MG  (/,  >h',  fh*) 
end ; 

applying  the  muliigrid  algorithm 

procedui'e  MG  (/,  u,  g) 

if  I I  then  u  •  -  SOLVE  (I,  u,  g) 
el  se 
beg  in 

for  t  --  1  to  n,  [while  ...]  do  u  *-  RELAX  (/,  u,  g); 

>  -  Ii  ^i-i  u; 

d  A i(g  -  Ah,u): 

for  i  •  -  I  to  n2  [while  . .  .  ]  do  MG  (1  -  1 ,  v,  d) ; 

u  •-  u  h  I,  -i «i(v  -  u); 

for  i  «—  1  to  n3  do  [while  ...]  u  <-  RELAX  (/,  u,  g) 
end ; 

After  ri !  relaxation  iterations  (procedure  RELAX)  have  been  performed  at  level  l,  MG  performs  a 
restriction  to  the  next  coaiscr  level  1-1.  It  then  calls  itself  recursively  on  the  coarser  level  n2 
times.  Finally,  it  performs  a  prolongation  from  the  co.u set  level  back  to  level  1.  following  up  with 
n3  more  iterations  or.  Ic*el  /.  I  lie  equations  on  the  coarsest  level  l  —  I  may  be  solved  to  desired 
accuracy  with  sufficiently  many  iterations  (procedure  SOLVE).  One  can  readily  show  that  w'hcn  MG 
is  invoked  on  icvel  X  it  uiu  RELAX  a  toL.i  of  n2x~'(.M  f  times  on  level  l  /  !  and  it  calls 
SOLVE  n2x_l  times  on  level  I.  In  general,  most  of  die  relaxation  iterations  arc  performed  on  the 
coarser  levels  |i Ictirker.  1980], 

The  optional  [while  ...]  clauses  denote  conditions  dial  may  be  checked  during  thc- 
computuiioii  and  used  to  terminate  some  iterations.  Dynamic  conditions,  oguvany  con.cigcucc 
rates  measured  by  crroi  norms,  are  incorporated  into  adaptive  coordination  schemes,  whereas  fixed 
schemes  arc  controlled  only  by  die  constants  nt,  n2,  and  u3  [lirandL  1977).  Although  adaptive 
schemes  tend  to  be  more  cllicient  in  practice,  fixed  schemes  lend  themselves  better  to  theoKtical 
analysis  and.  moreover,  they  a;e  easier  to  implement  on  distributed  local-mteiconncct  architectures 
due.  in  part,  to  the  absence  ol  crior  norms  which  require  global  compuLitions. 

3.  The  Lightness  Problem 

The  lightness  of  a  suif.icc  ;s  the  perceptual  correlate  of  its  reflectance.  Irradiancc  at  a  point 
in  the  image  is  proportional  to  t lie  product  ol  the  illuminance  and  retied. .nee  at  the  corresponding 
point  on  the  surface.  I  he  lightncw  problem  is  to  compute  lightness  from  image  irradiancc,  without 
any  precise  knowledge  of  either  reflectance  or  illuminance. 
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3.1.  Arudyals 

I  he  vctmex  theory  of  lightness  and  color  proposed  by  I  .and  and  McCann  {1971]  is  based  on  Ihe 
observation  that  illuminance  and  relkctancc  patterns  differ  in  their  spatial  properties.  Illuminance 
changes  arc  usually  gradual  and.  thcrelorc.  typically  give  rise  to  smooth  illumination  gradients, 
while  reflectance  changes  tend  to  he  sharp,  since  they  often  originate  from  abrupt  pigmentation 
changes  and  surface  occlusions.  Horn  [1974]  proposed'  a  two-dimensional  generalization  of  the 
1  .and-McCann  algorithm  for  computing  lightness  in  Mondrian  scenes,  consisting  of  planar  areas 
divided  into  subregions  of  uniform  matte  reflectance. 

I  .ct  'f(x,y)  be  the  reflectance  of  the  surface  at  a  point  corresponding  to  die  image  point  [x,y) 
and  let  S’x,y)  be  die  illuminance  at  that  point.  The  irradiancc  at  the  image  point  is  given  by 
h'(x,  y)  —  S{s.,y)x  lt{x,  y).  Denoting  the  logarithms  of  die  above  functions  as  lowercase  quantities, 
we  have  «[x,  y)  --  t  r(x, y).  Applying  the  Lnplaciati  operator  A  gives  d{x,y)  si  Ar {x,y)  => 

A.h(x,  y)  +  Ar(x,  //).  In  a  Mondrian,  illuminance  is  assumed  to  vary  smoothly  so  that  Aa(x,  y)  is  finite 
everywhere,  while  Ar(x,  y)  exhibits  pulse  doublets  at  intensity  edges  separating  neighboring  regions. 
A  dircsholding  operator  T  can  be  applied  to  discard  the  illuminance  component:  7’[<J(x,  y)]  — 
Ar(x,  y)  -  ■  J[x,y).  Hence,  the  reflectance  It  is  given  by  the  inverse  logaridim  of  die  solution  to 
Poisson's  equation 

Ar(x,  y)  --  f(x,  y),  in  D, 
where  fl  is  die  planar  region  covered  by  the  image. 

Horn  solved  the  above  partial  differential  equation  by  convolution  with  the  appropriate 
Green’s  function.  Wo  instead  pursue  a  local,  iterative  solution  based  on  die  finite  difference 
mcdiod.  Suppose  that  11  is  covered  by  a  uniform  square  grid  with  spacing  h.  We  can 
approximate  Ar  -  -  r,x  i  ryv  using  the  order  h'1  approximations  rxz  —  (tf+liJ  -2^  i  r[Li, j)/h? 
and  r,lv  (rj*  •  ( ,  t  r* •..,)//«*  to  obtain  a  standard  discrete  version  of  Poisson's  equation 

(r?+i,i  H  x\-ij  1  l<‘.;  m  1  r!‘,j  i  •  ~  u.j-  Hiis  denotes  a  system  of  linear  equations  with 

sparse  coefficient  matrix. 


Rearranging,  the  Jacobi  relaxation  step  is  given  by 


("■>')  I 


Wh  «*»)  ,  rh  <'*),  .*.  *"1  _  /,Jf*  .'l 

,Tr.  n.y  +*»-«. j  hr«.j  +  :  +ri.i-i  n 


where 'the  bracketed  superscripts  denote  the  iteration  index.  Jacobi  relaxation  is  suited  to  parallel 
synchrotuns  hardware,  whereas  the  Gauss-Seidcl  relaxation  step  given  by 


fc  ("+')  I  fjt 


1  (rh  !’*)  .  ..I,  !"H)  x  <”)  .  ,.h  |nM,_A2rO 

1  '<  -i,j  Kon  h  'i.i) 

is  more  suitable  on  a  serial  computer  and,  moreover,  requires  less  storage. 
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Figure  3.  Synthesized  Mondrian  imai/ox.  Huso  images,  input  to  the  algorithm.  cotiUiin  patches  of  uniform 
retlccumcc  anil  a  letWo-nghl  illumination  gradient.  I  lie  three  smaller  images  are  increasingly  coarser  sampled 
versions  of  the  largest  image  winch  is  tint  tao  pixels,  quantized  to  256  irradiance  levels. 

Wc  note  in  passing  that  Poisson’s  equation  Ar  /  is  the  Kulcr*  Lagrange  equation  for  the 
variational  principle  associated  with  a  membrane  problem,  'lire  solution  can  be  characterized  as 
the  deflection  v(x,y)  --  r(x,y)  of  a  membrane  subject  to  a  load  f(x,y).  and  it  minimizes  the 
potential  energy  functional  5(1;)  //nJ(«*  •*  nj)  -  fvdxdy  [Cotirant  &  Hilbert,  1953],  Blake 

[1984]  oilers  an  alternative  variational  principle  for  lightness.  Posing  die  '.gluncss  problem  as 
a  variational  principle  permits  the  direct  application  of  the  f'.nte  element  discretization  method, 
which  for  instance  does  not  require  a  uniform  discretization  of  fl. 

3.2.  Results 

A  four  level  aailui evolution  lightness  algorithm  (with  grid  sizes  129  v  129,  fifi  x  05,  33x33,  and 
17  x  17)  was  tested  on  a  synthesized  Mondrian  scene  consisting  of  patches  of  uniform  reflectance, 
subjected  to  an  illumination  which  increases  quadrntically  from  left  to  tight,  Hie  original  image, 
which  is  129  ;<  129  pixels  in  size,  and  three  co.aser-sampled  versions  ate  shown  in  P'ig.  3.  All 
images  arc  quantized  in  256  irradiance  levels.  The  grid  function  fl*  ..  diown  in  l:ig.  4,  was 
computed  by  maintaining  only  die  peaks  in  the  l.aplacian  of  r(‘;  Zero  boundary  conditions  were 
provided  around  the  edges  of  the  images,  and  the  compulation  was  started  from  the  zero  initial 
approximation  rj*,  0. 


8 


,» '■  *£i*n*m  >  ’  »  'ill'll 


1!  ::/0'  o:  ! . 


Ml  I  ilUKlh  «l  I  AXVIION  Mi  llions 


rigtiri  I.  1  •  •  . 

trr'tlic  I  apt  h.m  • 

l-'iti.  5  • 

RcCOnsiMi.t  . 

mi.il  : 1 1 s 1 1 » i ■  *:  >f 

62.  and  in.  k 

compute  .  • 
require-.  !. 
infor :•  i  i* -i  ■ 
lightness  i*:-  . 
coarser  v  ic 


;.i.  i  level.  i  !h-  o  functions  wire  obtained  by  maintaining  only  the  peaks 

V;  . . . 

.•Mi!  \*<>!:,ii i.*it  which  now  lacks  most  of  the  illumination  gradient. 

•  the  too, tions  vh,.wn  in  Fig.  4  required  .13.97  work  units.  The 

•  o:i  c.vh  leu'!  I'toiu  coarsest  to  finest  respectively  is  M2,  100, 
•  •I.  a  lose!  lightness  algorithm  required  about  500  work  units  to 

\  at  the  tine -t  level  in  isolation.  The  single-level  algorithm 

•  iOo-.n  i-«r  convergence  ,.s  there  are  nodes  across  the  surface,  since 

oidv  in  its  vnv-t  neighbors  in  one  iteration.  The  multilevel 
.  <!,.  at  k.  ..a  ,e  i!  I'iiip.iiMtcs  information  more  effectively  at  the 


4.  The  S  '•  k ; a <i  n ci  Problem 

In  ?:•  me  ,  •.•(»»  v.  m-  •aai.r..-.,  geometry,  scene  illuminance,  surface 

relleetancc  ..  •'  -  !!:••  Tape  tuun  shading  problem  is  to  recover  the  shape  of 

surfaces  I'r  »;:•  e,  -  it*  ..  .sHi-m;;  that  illuminance.  reflectance,  and  imaging  geometry 

arc  const.!  .t  uia..  .eve  ran  be  related  directly  m  surface  orientation. 

4.1.  Annlyj.m 

l.ct  i  i.r  :> :  t-  •  -  it. '  v  :;ii  '.< o' -t.uit  albedo  defined  over  a  bounded  planar  region  ft. 
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Figure  5.  The  lecnnsn  acted  M  >i:dna:i.  I  his  is  the  solution  computed  after  33.97  work  units  by  the  four-level 
lightness  algorithm.  Most  of  the  illumination  gradient  tn  i  tg.  3  has  been  eliminated. 

1  he  relationship  between  the  surface  orientation  at  a  point  (x,y)  and  the  image  irradiancc  tlicre 
E[i,  y)  is  denoieJ  by  li{p,  y).  where  p  -  Ux  and  </  -  are  the  first  partial  derivatives  of  the  surface 
function  at  (^.y)  The  shape-from-shading  problem  can  he  posed  as  a  nonlinear,  first-order  partial 
differential  equation  in  two  unknowns,  called  die  image-irradiaticc  equation:  h'(x,  y)  -  R(p,q)  =  0 
[Horn,  1975).  Surface  orientation  cannot  be  computed  strictly  locally  because  image  irradiancc 
provides  a  single  measurement,  while  surface  orientation  has  two  independent  components.  The 
image  irradiancc  equation  provides  only  one  explicit  constraint  on  surface  oncntalion. 

Iketichi  and  Horn  [1981]  proposed  an  additional  surface  smoothness  constraint  and  the 
use  of  surface  occluding  commits  as  boundary  conditions.  Since  die  ;>-< ;  parameterization  of 
sui face  orientation  becomes  unbounded  at  occluding  contours,  however,  surface  orientation  was 
rcparamcicii/ed  in  icons  of  the  (horn. Jed)  stenographic  mapping:  /  —  2 ap,  g  =  2 ay.  where 
a  -  1/(1  -  v  1  1  V1  <]')■ 

These  considerations  arc  fornt.ili/cd  by  a  v.ui.itional  principle  involving  the  minimization  of 
the  functional 

€(1,9)  -  /  J^ifl  H  ft)  4)  <!x  (ly  +  2  f  JjE{x,g)  At/,  g)\ltlr.dy. 

The  fust  integral  incorporates  the  surface  smoothness  constraint.  1  he  second  is  a  least-squares  term 
which  coerces  die  solution  into  satisfying  the  image  triad:. mcc  equation  by  Healing  the  equation  as 
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Figure  fi.  I  m  -plwic  images.  Ihcse  synthetic  images  input  to  the  algorithm  show  a  lainbeitian 

spline  distantly  ii:'!iiiiieii:'ii  thmi  die  viewing  direction.  the  three  smaller  images  are  increasingly  coarser 
sampled  vei>i>ms  >  i  the  I.iii  cm  image  which  is  120  >.  tetl  pixels,  quantized  to  ?5h  irradiance  levels. 

a  penalty  constraint  weighted  by  a  factor  X.  Other  variational  -formulations  for  shape-from-shnding 
have  been  suggested,  e.g...  [It  rooks  &  Horn,  1984J. 

The  I'tilei-  i  iriatiite  equations  are  given  by  the  system  of  coupled  partial  differentia!  equations 


A/  -  Xf/'’(x, y)  !((/,,, )\Rf  -  0, 

-  x;/s'(r.»)  -  «(/.»)]«*  *=o. 

Disereti/inc  .  iMH.uiims  on  a  uniform  grid  with  spacing  l>  using  standard  finite  difference 
approximation-.  ••!  F  the  Jacobi  relaxation  scheme 


where  ‘J'X'.j. 


t’»  »  I 


'•«  .  ; 


•%;*.,] '’0  1  *!'v,  .v(C./’“.^,”')j!w,JS). 


are  toe.  I  uw-t.'j  ■ 
<>!(/;>/.  and  !>'. 
relaxati  m  a:  ••• 
metnoi  reqmre! 
the  image. 


;  'l-l.;  '  1 1.  J  -  I  ’  -UHI  ,/‘g1‘>J]  S  g,  ,  IJ  t  P,,y  1  Ri,;'  f  ll/’l 

/'*  and  .,*■  at  node  (i,j)  (a  factor  of  1/4  has  been  absorbed  into  X),  11/  -■= 


On  a  sequential  computer,  we  prefer  to  use  the  analogous  Gnuss-Seidcl 
-  .ihile-.d  algorithm,  due  to  its  greater  stability,  faster  convergence,  and  reduced 
■nr--  Ai-piopri ate  boundary  comlitions  can  be  specified  at  occluding  contours  in 
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ITic  solution  computed  at  the  four  levels  after  6.125  work  units  arc  shown  in  Pig.  7.  'lire 
total  number  of  iterations  performed  on  each  level  from  coarsest  to  finest  respectively  is  32,  10,  4, 
and  4.  In  comparison,  a  single-level  algorithm  required  dose  to  200  work  units  to  obtain  a  solution 
of  llic  same  accuracy  at  the  finest  level  in  isolation.  As  in  die  ease  of  the  lightness  problem,  the 
single-level  algorithm  requires  at  least  as  many  iterations  for  convergence  as  there  arc  nodes  across 
the  surface,  since  information  at  a  node  propagates  only  to  its  nearest  neighbors  after  each  iteration. 
Convergence  is  somewhat  faster,  however,  because  shading  information  is  available  at  every  node 
inside  die  occluding  contour  to  constrain  surface  shape  according  to  die  image  irradiance  equation. 
In  any  case,  the  multilevel  shape-from  shading  algorithm  is  again  much  more  efficient  because  it 
enables  information  to  propagate  quickly  at  die  coarser  scales. 

To  obtain  a  representation  of  the  surface  in  depth,  die  surface  normals  in  Fig.  7  were 
introduced  as  orientauon  constraints  to  a  four-level  surface  reconstruction  algorithm  with  identical 
grid  sizes  (Tcrzopoulos,  1984a],  The  normal  vectors  were  first  transformed  from  the  f-g 
stenographic  parameterization  used  in  die  shape-from-shading  algorithm  to  die  p-q  gradient  space 
parameterization  used  in  the  surface  reconstruction  algorithm  using  die  formulas  p  =  -4 / /(/ 2  + 
t/2  -  4)  and  q  =  -4 g/{f2  -  g'“  -  4).  Nodes  outside  the  occluding  contour  of  die  sphere  were  treated 
as  depth  discontinuities.  Fig.  8  shows  die  surfaces  generated  by  the  algorithm  at  the  three  coarsest 
resolutions.  Hie  reconstruction  required  an  addiuonal  8.8  work  units. 

5.  The  Optica!  Flow  Problem 

Optical  flow  is  die  distribution  of  apparent  velocities  of  irradiance  patterns  in  die  dynamic 
image.  The  velocity  field  and  its  discontinuities  can  be  an  important  source  of  infomiation  about 
the  configurations  and  motions  of  visible  surfaces,  flic  optical  flow  problem  is  to  compute  a 
velocity  field  from  a  temporal  series  of  images. 

5.1.  Analysis 

Horn  and  Schunck  ]  19S 1  ]  suggested  a  technique  for  determining  optical  flow  in  the  restricted 
ease  where  the  observed  velocity  of  image  irradiance  patterns  can  be  attributed  directly  to  small 
interframe  motions  of  surfaces  in  the  scene.  Under  these  circumstances,  the  change  in  image 
irradiance  at  a  point  (c,  </)  in  the  image  plane  at  time  t  and  die  motion  of  the  i r radiance  pattern 
can  be  icluied  by  die  flow  equation  Ltu  -r  Ky v  i  Kt  --  0,  where  IC(r,y,t)  is  die  image  irradiance, 
and  u  -  dx/iU  and  v  dy/di  arc  the  optical  flow  component  functions. 

An  additional  constraint  is  needed  to  solve  this  linear  equation  for  the  two  unknowns  u  and  v. 
If  opaque  objects  undcigo  rigid  motion  or  deformation,  most  points  have  a  velocity  similar  to  that 
of  their  neighbors,  except  where  surfaces  occlude  one  another.  Observing  that  the  velocity  field 
varies  smoodily  almost  everywhere,  optical  flow  can  he  determined  by  finding  the  flow  functions 
u(x,y)  and  v(t,  y)  which  minimize  the  functional 
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C(u,v)  =  a 2  J  J  ( u]  -  u‘)  +  (v*  -t-  vl)dzdy  +  J  f  V'*u  +  Evv  +  Etf  dxdy, 

where  a  is  a  constant.  The  first  term  is  the  smoothness  constraint,  while  the  second  is  a  least-squares 
penalty  expression  which  coerces  the  How  field  into  satisfying  the  flow  equation.  Related  variational 
formulations  of  the  optical  flow  problem  have  been  suggested  (e.g.,  [Nagel,  1983],  [Cornelius  and 
Kanadc,  1983]). 

The  F.ulcr-l.agrange  equations  for  the  functional  £  arc  given  by  [Horn  and  Sch uncle.  1981] 

E\u  +  EzEyv  --  «sAti  -  EzEt, 

Ez  Eyii  -r  Eyii  — ■  ft  ^  A v  —  EyEf 

Assuming  a  cubical  network  of  nodes  with  spacing  h ,  where  t.  j.  and  k  index  nodes  along  ihc  x,  y, 
and  t  axes  respectively,  v.c  use  the-  following  finite  difference  fonnulas  to  discreti/.c  the  differential 
operators 


;  /  ’  ]h  _ ®  /  t/K  j,»X  V 

I* 21-  v'i  *  I  "  *'i  -  I 

=  2/i + 

Ahu  =--  ^  (h;u,hiJ  Jt]  -  u]'j^)| 

>:^,1  ~  v*,,). 


where  0[u,hi>it]  =■■  m.*  J  i.*  +  u?+IiJi*  J  and 

0[vj*,  J  =_  .{(V*.,  ,  .  r  v|*  , ih  i  vl\,  ,  *  -r  v]'J.  ,  k).  Oilier  approximations  arc  possible,  including 
those  suggested  by  Horn  and  Sclumck  which,  however,  require  over  four  limes  the  computation  per 
iteration  to  gain  some  improved  attenuation  of  high  frequency  error.  Given  dynamic  images  over 
at  Rest  three  frames,  a  symmetric  central  difference  formula  -  E*}  k_ ,) 

would  be  prefeiablc.  provided  it  is  stable. 

Substituting  flic  above  approximations  into  the  lailer- 1  agrangc  equations  and  solving  for 
u[\  k  and  vj1  t  yields  the  following  Jacobi  relaxation  formulas 
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where  /«*,.*  -  ([E S  ■  (  /'T ^  * )2  h-  ~  <r 

,  --  [/■V,h;y/vu,hJ  J i'-Al;. yP\\*j a  ■  .  riienaiur.il  boundary  conditions  of  tlic  zero 

normal  derivative  arc  appropriate  on  ihc  boundaries  of  surfaces,  t  hey  can  be  enforced  by  copying 


values  to  boundary  nodes  from  neighboring  interior  nodes. 
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Figure  9  Lambertian  sphere  images.  These  synthetic  images  input  to  the  algorithm  at  four  resolutions  depict 
a  uniformly  expanding  lambcnian  sphere,  disumtly  illuminated  from  the  viewing  dircoion.  Frames  for  the 
first  time  instant  are  shown  to  the  left  of  frames  for  the  second  tune  instant. 

5.2.  Results 

A  four  level  optical  flow  algorithm  (with  grid  si/.cs  129  X  129.  65  x  65,  33  X  33,  and  17  X  17) 
was  tested  on  a  synthetically-generated  image  of  a  Lambertian  sphere  distantly  illuminated  from 
the  viewing  direction  by  a  point  source.  The  sphere  expanded  uniformly  over  two  frames.  The 
first  frame,  which  is  129  x  129  pixels  in  si/e,  and  three  toarscr-samplcd  versions  arc  shown  in  the 
left  half  of  Fig.  9.  T  he  next  frame,  in  which  the  sphere  has  expanded  is  shown  in  the  right  half 
of  the  figure  All  images  arc  quantized  to  256  irradiance  levels.  T  he  velocity  field  was  specified 
around  die  occluding  contour  of  the  sphere,  and  by  treating  the  contour  as  a  possible  How*  field 
discontinuity,  «  and  v  were  allowed  to  make  discontinuous  transitions  across  it.  T  he  computation 
was  started  from  die  zero  initial  approximation  u  —  v  —  0. 

The  solution  computed  on  the  three  coarsest  levels  after  4.938  work  units  arc  shown  in  Fig. 
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Figure  10.  Velocity  vectors  for  ihc  expanding  UiinFeriiaii  sphere.  Hie  solution  at  the  three  coarsest  resolutions 
thui  were  obtained  after  4.93h  woik  unite  arc  shown  (the  finoM-lcvcl  solution  is  loo  <len»e  to  depict).  _ 

10  as  velocity  vectors  in  ji/'spacc.  I  lie  total  fine. her  ol  iterations  pci  tunned  on  each  level  from 
coarsest  to  finest  respectively  is  40,  5,  4,  am'  .V  In  comparison,  a  single-level  algorithm  required 
37  woik  units  to  obtain  a  solution  of  the  same  accuracy  at  the  finest  level  in  isolation.  Again, 
the  multilevel  algorithm  is  more  cilicicm  because  it  propagates  information  quickly  at  the  co<*rscr 
scales.  Ole /ci  [1934]  also  reports  improvements  consistent  with  ours  with  regard  to  the  convergence 
rate  of  a  multilevel  optical  How  algorithm  relative  to  a  single  level  algorithm,  lie  employed  the 
I  lom-Sclu  nek  iclaxalion  formulas  for  Ins  implementation. 
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6.  Multigrid  Methods,  Regularization,  and  Stochastic  Relaxation 

A  primary  purpose  of  low-level  visual  processing  is  to  reconstruct  relevant  physical  charac¬ 
teristics  of  3-D  scenes  from  their  images.  Wc  have  considered  in  this  paper  three  different  visual 
reconstruction  problems  —  the  compulation  of  lightness  from  an  image  (a  2-D,  static  reconstruction 
problem),  shape-from-shading  (a  3-D,  static  problem),  and  optical  flow  (a  2-D,  dynamic  problem). 
It  was  possible  to  apply  multigrid  methods  because  each  of  these  problems  was  formulated  as  a 
variational  principle  or  associated  partial  differential  equation. 

As  inverse  mathematical  problems,  visual  reconstruction  problems  tend  to  be  mathcinatical'y 
ill-posed,  in  that  existence,  uniqueness,  and  stability  of  their  solutions  cannot  be  guaranteed  a 
prion  fPoggio  and  Torre,  1984],  Among  the  systematic  techniques  that  have  been  developed  to 
tackle  ill-posed  problems  is  the  method  of  regularization  [Tikhonov  and  Arsenin,  1977).  Ihrough 
regularization  analysis,  ill-posed  visual  problems  can  be  restated  as  well-posed  variational  principles 
by  restricting  the  possible  solutions  with  appropriate  stabilizing  functionals.  In  general,  the 
smoothness  properties  of  stabilizers  must  be  controlled  near  discontinuities  [Terzopoulos,  1984b). 
1  .  fingly,  the  same  stabilizer  was  used  to  impose  the  smoothness  constraint  in  both  the  shape- 

from-shading  and  optical  flow  problems. 

A  major  attraction  of  regularization  analysis  is  that  it  leads  systematically  to  variational 
principles  which  permit  advantageous  use  of  multigrid  relaxation  methods.  As  a  visual  algorithm 
design  strategy,  regularization  analysis  applied  in  conjunction  wdth  multigi  id  methoc ology  promises 
to  impact  on  a  bioadcr  spectrum  of  visual  reconstruction  problems,  including  image  reconstruction 
and  discontinuity  detection  [Gcman  and  Goman,  1984),  stcrcopsis  [Marr  and  Poggio,  1977), 
registration  [llajesy  &  Broil,  1982],  motion  field  interpolation  [Hildreth,  1984],  shape  from-contour 
[Brady  &  Yuillc,  1933],  and  siruclurc-from-motion  [Ullman,  1979). 

An  issue  of  concern  is  that  die  regularization  of  visual  reconstruction  problems  cannot  always 
be  expected  to  lead  to  convex  variational  principles  having  a  unique  absolute  extremum,  without 
relative  extrema.  Unfortunately,  classical  relaxation  or  gradient  descent  methods  are  not  directly 
applicable  to  nonconvcx  variational  principles,  since  they  often  get  trapped  in  relative  extrema 
Stochastic  relaxation  algorithms  (such  as  simulated  annealing)  do  not  suffer  this  disadvantage 
[Kirkpatrick  cl  al,  1983;  Hinton  &  Scjnowski,  1983).  Nonetheless,  since  stochastic  relaxation 
searches  for  absolute  extrema  with  processors  Uiat  arc  restricted  to  local  interactions,  it  too  suffers 
serious  inefficiencies  in  propagating  constraints.  The  inherently  slow  convergence  rates  arc  further 
aggravated  by  the  nondctcrministic  nature  of  die  local  computations.  Multigrid  methods  may 
ameliorate  these  problems  by  facilitating  constraint  propagation  ihrough  the  use  of  coarser  scales. 


7.  Conclusion 


Many  important  problems  in  low-level  computer  vision  can  be  formulated  as  variational 
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principles  or  as  partial  differential  equations.  A  particular  source  of  such  formulations  is  the 
regularization  analysis  of  til-posed  visual  reconstruction  problems.  Once  discretized,  variational  and 
differential  formulations  arc  amenable  to  numerical  solution  by  iterative  relaxation  methods,  which 
readily  map  into  massively  parallel  computer  architectures.  However,  distributed  local-support 
computations  arc  inherently  inefficient  at  propagating  constraints  over  the  large  network  or  grid 
representations  that  arc  encountered  in  computer  vision  applications. 

In  our  previous  work  on  surface  reconstruction  algorithms,  we  established  that  multiresolution 
relaxation  techniques  can  overcome  this  inefficiency,  without  sacrificing  the  local-interconnect  nature 
of  the  coi  iputations.  T  his  has  been  corroborated  in  the  present  paper  by  successfully  applying 
muliigrid  methods  to  die  well-known  problems  of  computing  lightness,  shape-from-shading,  and 
optical  flo a  from  images.  The  novel  nu.liircsolution  algorithms  that  we  designed  in  the  context 
of  each  of  these  problems  were  show  n  to  be  substantially  more  efficient  than  the  published  single 
level  versions. 

Beyond  its  effectiveness  as  a  (local)  convergence  acceleration  strategy,  our  adaptation  of 
muliigrid  ncihouology  also  leads  to  iterative  algorithms  that  compute  mutually  consistent  visual 
rcprcsenta.ions  over  a  range  of  spa lia!  scales.  Multiresolution  representations  appear  to  be  crucial 
in  interfacing  low-level  visual  processing  to  subsequent  tasks  such  as  recognition,  manipulation, 
and  navigation. 
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